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Abstract 

The structure of a genetic network is uncovered by studying its response to ex- 
ternal stimuli (input signals). We present a theory of propagation of an input signal 
through a linear stochastic genetic network. It is found that there are important 
advantages in using oscillatory signals over step or impulse signals, and that the 
system may enter into a pure fluctuation resonance for a specific input frequency. 

The nature of a physical system is revealed through its response to external stimulation. The 
stimulus is imposed upon the system and its effects are then measured, Fig. 1(a). This approach 
is widely used in biology: a cell culture perturbed with a growth factor, a heat shock etc. The 
data measured contain the initial information encoded into the stimulus plus the information about 
the intrinsic characteristics of the system. The more parameters the experimentalist can adjust 
to craft the perturbation stimulus, the more information about the system can be revealed. In 
recent years we witnessed a tremendous increase in measurement capabilities (e.g., microarray and 
proteomic technologies, better reporter genes). However, the success of the systems approach to 
molecular biology depends not only on the measurement instruments, but also on an effective design 
and implementation of the input stimulus, which has not been thoroughly explored. Traditionally, 
two types of time dependent stimuli are at work in molecular biological experiments pQ, For 
example a step stimulus is obtained when at one instant of time a growth factor is added to the 
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medium, graph (a) in Fig. 1(a). The stimulus from Fig. 1(a) graph (b), is a superposition of two 
step stimuli. The investigator can control the height of the step stimulus (the concentration of the 
growth factor) or the time extension of the heat shock. The cells respond to these stimuli only 
transiently. The response is dampened after some time and becomes harder to detect it from noise. 
To overcome the noise, the concentration of the stimulus is typically increased to the point where 
the strength of the stimulus raises far above its physiological range. 

We propose to implement a molecular switch at the level of gene promoter and use it to impose 
an oscillatory stimulus. In the absence of experimental noise, any stimulus can be used to determine 
the genetic network input-output properties. However, in the presence of experimental noise, the 
oscillatory input has many advantages: (1) the measurements can be extended to encompass many 
periods so the signal-to-noise ratio can be dramatically improved; (2) the measurement can start 
after transient effects subside, so that the data becomes easier to be incorporated into a coherent 
physical model; (3) an oscillatory stimulus has more parameters (period, intensity, slopes of the 
increasing and decreasing regimes of the stimulus) than a step stimulus. As a consequence, the 
measured response will contain much more quantitative information. Experimental results from 
neuroscience prove that oscillatory stimulus can modulate the mRNA expression level of genes. For 
example, c-fos transcription level in cultured neurons is enhanced 400% by an electrical stimulus 
at 2.5 Hz and reduced by 50%at 0.01 Hz, Also, the mRNA levels of cell recognition molecule 
LI in cultured mouse dorsal root ganglion neurons change if the frequency of the electric pulses 
is varied. The expression level of LI decreases significantly after 5 days of 0.1 Hz stimulation but 
not after 5 days of 1 Hz stimulation [I]. To extend the oscillatory approach to other type of cells, 
a two-hybrid assay, [Sj, can be used to implement a molecular periodic signal generator, Fig. 1(c). 
The light-switch is based on a molecule (phytochrome in [2j) that is synthesized in darkness in the 
Ql form. When Ql form absorbs a red light photon ( wavelength 664 nm) it is transformed into 
the form Q2. When Q2 absorbs a far red light ( wavelength 748 nm) the molecule Q goes back 
to its original form, Ql. These transitions take milliseconds. The protein P interacts only with 
the Q2 form, recruiting thus the activation domain (AD) to the target promoter. In this position, 
the promoter is open and the gene is transcribed. After the desired time elapsed, the gene can be 
turned off by a photon from a far red light source. Using a sequence of red and far red light pulses 
the molecular switch can be periodically opened and closed. 
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Figure 1: (a): The genetic network response depends on the type of the applied stimulus, (b): 
An autoregulatory network. The gene G is under the influence of a cofactor C that rhythmically 
modulates the activity of the promoter P. The matrix H contains the parameters that dictates the 
transition probabilities of the stochastic model. The transition probability per unit time from r to 
r + 1 mRNA molecules, T(D,r,p; D,r + l,p;t),gs modulated by the oscillatory signal generator. 
The DNA, D, and the protein, p, do not change in this transition, (c), adapted from [5]. The gene 
is turned on with a red light pulse of wavelength A = 664 nm. With a far red light of wavelength 



There are four input parameters that can be varied: the period (T), the time separation between 
the pulses (s), and the amplitude (A) of the red and the far red pulses. The mRNA concentration 
profile will depend on these parameters and can be measured with a high throughput technology UJ. 
Protein levels will also depend on the input signal. The proteins can be recorded with 2D PAGE 
analysis or mass spectrometry. If one single gene product is targeted, than a real-time luminescence 
recording can be employed [7]. A periodic generator can be used to investigate biological system 
for which the mRNA and protein concentrations naturally oscillate in time. An example of such 
a system is the circadian clock that drives a 24-hour rhythm in living organisms from human to 
cyanobacteria. The core oscillator is a molecular machinery based on an autoregulatory feedback 
loop involving a set of key genes (Bmal, Perl, 2, 3, Cryl,2 , etc.) jSj. Experimental procedures used 
to elucidate the clock mechanism are based on measuring the circadian wheel-running behavior of 
mice under normal light/dark (LD) cycles or in constant darkness (dark/dark or DD) conditions. 
Experimental evidence demonstrates that laws of quantitative nature govern the molecular clock. 
For example, [§], the internal clock of cry 1 mutants have a free-running (i.e. DD conditions) 
period of 22.51 ± 0.06 h which is significantly lower than the period of a wild-type mice which is 
23.77 ± 0.07 h. Quite opposite, a cry 2 mutant have a significantly higher period of 24.63 ± 0.06 h. 
In LD conditions, both mutants follow the 24 h period of the entrained light cycles. A double cryl,2 
mutant is arrhythmic in DD conditions and follow a 24 h rhythm in LD conditions. To explain 
these experimental values we suggest using a light switchable generator to drive the expression 
level of cryl, 2 and measure the dynamics of transcription and the translation for the rest of the 
key clock genes. Another application of the periodic generator is to modulate a constitutively 
expressed gene by superimposing an oscillatory profile on top of its flat level. Then, the genes 
that show a modulation with a frequency equal to the generator's frequency will be detected by 
a microarray experiment. Why is this approach different from the one where a step stimulus is 
used? Because the frequency of the generator is not an internal parameter of the biological system. 
The genes that interact with the driven gene will be modulated by the input frequency. The 
rest of the genes will have different expression profiles, dictated by the internal parameters of the 
biological system. This point of view is supported by our findings, that the circadian clock 
(which is an endogenous periodic signal generator) propagates its output to only 8 — 10% of the 
transcriptome in mice peripheral tissues (liver or heart). In contrast to the oscillatory input, when 
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a step stimulus is applied, all the expression profiles are dictated by the internal parameters of the 
biological system. Except for the height of the step stimulus (the dose of the factor applied) there 
is no external parameter implemented into the input signal. As such, it is difficult to separate 
those genes that directly respond to the input signal and to consequently avoid artifacts. With the 
applications described in mind, we study the propagation of an input signal through a stochastic 
genetic network. 

1 The response of a stochastic genetic network to an input stim- 
ulus 

The effects of an oscillatory input were previously studied on specific examples using models based 
on differential equations [E| , [15 j [121 ■ The stochastic character is embedded into these equations 
as an exterior additive term. In contrast, we compute the generator's effects on the mean and 
fluctuation of the gene products using a stochastic model pQ, [H|,[2]- in this way, the generated 
stimulus and the noisy nature of the cell are entangled in the stochastic genetic model. For a 
network of n genes the state of a cell is described by the mRNA and protein molecule numbers: 
q = (ri, ...,r n ,pi, ...,p n ). We assume that, during any small time interval At, the probability for 
the production of a molecule of the ith type is (X)|=i AijQj + Gi(t))At, i.e. qi is increased by 1 with 
the above probability. The function Gj(t) represents the time varying input signal and modulates 
the mRNA production only: G = (g\(t), . . . , g n (t),0, . . . , 0) T (the superscript T is the transposition 
operation that transforms G into a column vector for notational convenience in what follows). The 
parameter Ay represents the influence of the jth type of molecule on the production rate of a 
molecule of the ith type. Similarly, there is a matrix of parameters governing the degradation 
rates of the molecules. For simplicity, we assume that the input stimulus directly affects only the 
production rates. The mean \i = (q) and the covariance matrix v = ((q)) = ((q — (q)){q — (l}) T ) 
of the state q are driven by the generator G. 

The transfer of the signal from the generator through the genetic network to the output mea- 
sured data is encapsulated in a set of transfer matrices. Specifically, let H = A — V and denote 
the Laplace transforms of fi and G by L\i and CG. Here and in what follows, fi and G are repre- 
sented as column vectors. The connection between the mean and the generators is given by formula 
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Figure 2: Response of a stochastic genetic network to an oscillatory input. The Laplace transform 
C change the dynamic variable from time to frequency. In the vec(X) all the elements of the matrix 
X are arranged in a column vector. 

(1) which is typical for a deterministic linear system. However, the genetic system is stochastic 
and the measure of the intrinsic noise is quantified by the covariance matrix u. The effect of the 
stimulus generators is most transparent if we split v in a Poisson and a non-Poisson component: 
v = diag(fx) + X. Here diag(fj,) represents a matrix with the components of the vector /i on its 
diagonal, all the other terms being zero. For a Poisson process, X = and thus the term diag(ii) is 
called the Poisson component of v. The non-Poisson component X = v — diag{(j) can be expressed 
in terms of the generators (Appendix and Supplementary Material): 

1 



Cvec(X) = * [{1 ® H + H ® l)L + 2Lr]-!—£G . 

s — 1® H — H ®1 s — H 



(1) 
(2) 



The vec(X) is a vector constructed from the matrix X by stripping the columns of X one by 
one and stack them one on top of each other in vec(X). We emphasize here that the time variation 
of the generators G in (2) can take any form and is not bounded to be periodic or a step stimulus. 
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There are 3 matrices that transfer the information from the generators to the non-Poisson 
component, Cvec(X) = M 3 M 2 MiCG. The first, M 1 = (s - H^ 1 , is the same as the transfer 
matrix for the mean. The second, M 2 = (1 <8>H + H <g> 1)L + 2LT, breaks the symmetry between the 
degradation and production parameters that are otherwise hidden in the matrix H = A — F. The 
<£> is the Kronecker product of two matrices. The matrix L (with elements and 1) is the lifting 
matrix from dimension of the mean (2n values) to the dimension of vec(X) ( An 2 values). The 
third matrix is Ms = (s — 1® H — H ® 1) . If Aj are the eigenvalues of H then all combinations 
Aj + Xj are the eigenvalues of 1 <g> H + H <g> 1. Thus M3 represents the analog of Mi in the space of 
covariance variables. 

For a step stimulus, these eigenvalues are of primal importance: the measured signal is a super- 
position of components with different eigenvalues and has a complicated mathematical expression. 
However, for a periodic stimulus, the frequency of the external generator is the important parame- 
ter. This frequency is fixed by the experimentalist not by the biological system. Only the phase and 
the amplitude of the output signal depends on the system's eigenvalues and the mathematical form 
is less cumbersome then for the step stimulus. The input-output relations, (1) and (2), were derived 
from the Master Equation written for the probability of the states of the genetic network. Thus we 
must specify the initial conditions for the probability of the states. These conditions refer here to 
states for which one molecular component vanishes (<& = 0, for one i). The input-output relations, 

(1) and (2), are independent of these boundary states if the T matrix is diagonal. A diagonal T 
matrix was used in and we will use it also in the example that follows. Tools developed in the 
field of System Identification can be used to create models for the networks under study, [T^]. The 
difference between the System Identification classical models and a genetic network is that the later 
is a stochastic process by nature, whereas the former are deterministic models with a superimposed 
noise from external sources. However the formulas that describe the relations between the mean 
and covariance of the stochastic process and the input signals, (1) and (2), are of the same general 
nature as those used in System Identification Theory, [T2|. In the next section we will use (1) and 

(2) to analyze one of the most fundamental regulatory motifs in a genetic network: an autoregula- 
tory gene that acts upon itself through a negative feedback, |1).|18|.[T5]. The fluctuation can drive 
this biological system out of its equilibrium state. [2*0]. 
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2 Fluctuation resonance 

Four parameters characterize the system: the feedback strength A12 = —h ; the translation rate 
A21 = k p , and two degradation rates, F±i = j r , T22 = 7 P - The gene regulation is under the control 
of its own protein product and the protein activity is modulated by a cofactor. The cofactor is 
driven by a periodic light switchable generator g(t) = k^ + acos{uit), Fig. 1(a). Before the generator 
is applied, the transcription rate is equal with ko and the system is in a steady state. Through 
the transfer matrices, (1) and (2), the light generator will impose a periodic evolution of the mean 
and covariance matrix for mRNA and protein product. We denote the mean mRNA by {r(t)) 
and the mean number of protein by (p(t)). We will concentrate on the protein number in what 
follows. After the transients are gone, (p{t)) = P + Pie iult + Pfe~ iujt , that is the protein number 
will oscillate with an amplitude P\ on top of a baseline Poi here * represents complex conjugation. 
The fluctuation of the protein number, {(p(t))) , differs from the mean number by a quantity that 
we denoted by X pp (t): {{p(t))) = (p(t)) + X pp (t). For a pure Poisson process, {{p(t))) = {p(t)). 
Thus the term X pp (t) represents the deviation from a Poisson process. If there is some information 
about the genetic system that can be uncovered by measuring not only the mean but also the 
covariance matrix, then this information is hidden only in the non-Poisson component X pp (t). The 
quantity X pp (t) is not interesting only from a statistical point of view but also from a dynamical 
one. The equation for the time evolution of ((p(t))) takes its most simple form if it is written 
for X pp (t). That is, the time dependence of the mean value must be subtracted from the time 
evolution of ({p(t))) . Similar to the mean value, the non-Poisson component of the fluctuation 
will oscillate in time, X pp {t) = X p fi + X Pt ie tujt + X* 1 e~ tujt with complex amplitude X Pi i. The 
relative strength of the fluctuation versus the mean value can be described using the Fano factor, 
C0 : {(p(t))) I {p(t)) = 1 + ^pp(*)/ • For oscillatory inputs, the response of the network is best 
described in frequency domain rather than in time. In frequency domain, as an analog of the Fano 
factor we consider the ratio of the amplitude of X pp {t) versus the amplitude of {p(t)) . 



Here uj\ = 7 r + 7 p . The complex amplitudes X P) \ and Pi depend on the input frequency and 
therefore resonance phenomena can be detected in the system. If the light switchable generator 
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Figure 3: Fluctuation resonance. The amplitude X p> i of the non-Poisson component is much higher 
than the amplitude of the mean protein number, Pi, at u = 2ojq. 

oscillates with double the natural frequency uJq = hk p + 7 r 7 P , that is , u = 2u>o we find a state of 
resonance for fluctuation and not for the mean, Fig. 3. 

For co = 2ojq the system will be in a pure fluctuation resonance. In such a situation the molecular 
noise can drive the cell out of its equilibrium state, which can have dramatic consequence on the 
cell fate. Our model being linear cannot cover the entire phenomena that accompanies a system 
whose state is close to resonance. However, a linear model suggest the existence of pure fluctuation 
resonance. At fluctuation resonance, the deviation from a Poissonian process is high. The oscillation 
amplitude for protein fluctuation is much greater then the amplitude of the mean. Experimental 
results |21) show that typical values for the ratio k v j^ r are 40 for lacZ and 5 for lacA. This suggests 
that there are natural conditions for a strong height fluctuation resonance, Fig. 3. However, for a 
sharp fluctuation resonance (small half width), we need h > j r or j p , a condition that does not 
appear in all genetic networks. It is with the help of the experimental study that we will clarify 
why some biological systems can sustain fluctuation resonance and others not. Beside resonance, 
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the frequency response provides other insights into the structure of the autoregulatory system. The 
parameters of the system can be read out from the measured data. The frequency response of the 
mean values behave like the response of a classical linear system to input signals. The new aspects 
are those related to fluctuations. Like X pp {t) and X rr (t), the correlation coefficient between the 
mRNA and protein number will oscillate in time: X rp (t) = X rp> o + X Tp ^e lult + X* pl e~ lult with 
amplitude X rP) ±. Taking the ratios of the amplitudes: |X rPi i| 2 /|X rj i| 2 = (l/4h 2 )uj 2 + ^ 2 ./h 2 , 
\X rPi i\ 2 /\X P) x\ 2 = (l/4:k 2 )uj 2 + 7p/fc 2 , we observe that all four parameters of the system can be 
estimated from the slopes and the intercepts of the above ratios as a function of lo 2 . Detail formulas 
for each amplitude are given in the Supplementary Material. 

3 The spectrum, the experimental noise and the importance of 
the input stimulus 

We described the use of a periodic signal to decipher a genetic network. Traditionally, a step 
stimulus is employed in biology for pathway detection (i.e., adding a growth factor to the culture). 
From the response to a step stimulus we can extract, in principle, the parameters of the system. 
The natural question is then: why should we generate a periodic stimulus when there is already 
a step stimulus in use? Seeking an answer, we notice that the measured data in our studied 
example can be expressed as a sum of exponentially decaying functions, e~ xt , if a step stimulus 
was used (Supplementary Material). For a periodic input, the response contains only exponentials 
with imaginary argument, e tlJjt . Mathematically, the main difference between exponentials with real 
arguments, e~ xt , and those with imaginary arguments, e lujt , is that with the former we can not form 
an orthogonal basis of functions whereas such a basis can be formed with the later. If we depart from 
our example, we can say that in general, the response of the network to a step input will be a sum of 
components which are not orthogonal on each other. The time dependance of these nonorthogonal 
components can be more complex than an exponential function; they can contain polynomials in 
time or decaying oscillations, depending on the position in the complex plane of eigenvalues of the 
transfer matrix H. Contrary, the permanent response obtained from a periodic input is a sum of 
Fourier components that form an orthogonal set. Orthogonal components are much more easy to 
separate than nonorthogonal ones. This mathematical difference explains the advantage of using 
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oscillatory inputs. However, an argument can be made that increasing the number of replicates will 
be enough to recover the step response form noise. In what follows we study how many replicates 
we need to successfully fight the experimental noise. We will show that we need fewer replicates if 
the genetic network is probed with an oscillatory generator than with a step signal. To keep the 
argument simple, we will study the difficulty of separating nonorthogonal components for a network 
for which the response to a step stimulus is a sum of decaying exponentials. The argument can 
be extended to other types of nonorthogonal components, but this line of thought falls out of the 
scope of this paper. The measured data being a superposition of exponential terms can be written 
as: 

/(*)= Hs(x)K(xt)dx, (4) 

J Xl 

with K(xt) = e~ lxt for the periodic response and K(xt) = e~ xt for the step stimulus. The 
spectral function S(x) depends on the network's parameters and on the type of the input signal. 
For example, the spectrum of the autoregulatory system for a periodic input is S(x) = Sq5(x) + 
S\5(x — iuj) + S*S(x + iu), where 5(x) is the Dirac delta function. The coefficients So, S% take 
specific values if the spectrum refers to mean mRNA, proteins or their correlations. For example, 
for the protein fluctuation: 

o Y fc P 2 fcp (7 P - h) 7r 

Jo — ^ P ,o — 7 j W 

wo Wl 

S 1 = X i = — ^~ ilp + - + — — - 2 ilr ^ kp2 (6) 

Detailed description of the spectrum for an autoregulatory network is given in the Section 5 of 
the Supplementary Material. For oscillatory inputs that are not pure cosine function and for more 
complicated networks, the spectrum is more complex, but still is connected with the measured data 
like in (jlj). The spectrum S{x) carries information about the parameters of the genetic network 
and it can be recovered from the data f(t). The network's parameter can be estimated from the 
spectrum once a model of the network is chosen. Our goal is to show that the spectrum obtained 
from an oscillatory input signal is much less distorted by the experimental noise than the spectrum 
obtained from a step input. Laboratory measurements are samples of f(t) at N discrete time 
points. Given a finite number iV of measured data points, /i, ■ ■ ■ , /jv, the spectrum for the periodic 
case S(x) can only be approximated as a weighted sum of N terms, (Supplementary Material): 
S( x ) = J2k=i( s k + e k/ Pk)®k( x )- Each term, (s^ + e&/ ^)Ofc(x), contain a function @k( x ) that do 
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not depend on the measured data, and the weights Sfc + £fc//3fc that are computed from the measured 
data fx, ■ ■ ■ , /at. In the absence of experimental noise, e& = 0, all N coefficients Sk can be computed 
from the measured data. When experimental noise is present, 7^ 0, what we compute from 
measured data is Sk + ek/ /3k, an d we cannot separate Sk from it because we do not know the actual 
value for e^. The best we can do is to use only those terms for which Sk > tk/Pk, so the effect of 
the distortion on Sk is not large. Unfortunately, the distortion increases as f3k goes smaller, which 
actually happens when k increases. A term can be recovered from noise if (3k~ l < Sk/ek- Usually, 
this relation is valid for k = 1 • • • J p , with J p being the last term that can be recovered. A similar 
relation holds for the exponential case, with otk instead of /3k and J e instead of J p . It is desirable 
that both cutoffs (J p , J e ) be as close as possible to the number of sampled points, N. The striking 
difference between the two cases is that the cutoff J p is much larger then the cutoff J e . This is a 
consequence of the fact that the numbers ak decrease exponentially to 0, [S], whereas (3k stays close 
to 1 for many k before eventually dropping close to zero, [Sj. This huge difference between ak and 
[3k has its origin in the fact that the set of functions of time, exp(-Xt), indexed by A, do not form an 
orthogonal set, whereas the functions exp(iujt), indexed by u, are orthogonal. In theory, however, 
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Figure 4: How many replicates we need to recover a given spectral component. 

we can still hope that a step stimulus can deliver good estimates if the noise e& is reduced using r 
replicates (e& — > e^/y^r). This is not the case. Fig. 4 represents the number of replicates needed to 
recover the component J e or J p if the Signal to Noise Ratio is 10 (SNR = sj e /e j e = sj p /ej p = 10). 
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The number of replicates grows very fast in the exponential case (for SNR = 10 and N = 20 we 
need 269 replicates for the 4th spectral component), whereas in the periodic case, the number of 
replicates stays low for many spectral components (only for the 17th component it raises to 14, 
with SNR = 10 and N = 20). 

4 Conclusions and Discussions 

We studied the response of a linear stochastic genetic network to an input stimulus (signal). We 
provide a general formula that relates the mean and covariance matrix of mRNAs and proteins to the 
input generators. The particular type of periodic signals was studied in detail for an autoregulatory 
system. We found that fluctuation resonance can manifest in such systems. Besides interesting 
physical phenomena that can be detected using a periodic signal, the oscillatory input is useful 
for experimental noise rejection. We compared two experimental designs: one that uses a step 
stimulus as a perturbation and another one that uses a periodic input. We concluded that the 
response of the genetic network to a periodic stimulus is much easier to be detected from noise 
than the response of the same network to a step stimulus. This conclusion applies whenever the 
response of the network to an oscillatory input is a sum of Fourier components. This can be the 
case for many nonlinear networks. However, the input-output relations, (1) and (2), applies only 
to a linear stochastic model. A linear model is a good approximation around a steady state of 
the genetic network. A genetic network is a nonlinear system and can have several steady states. 
If the signal generator does not vary in time, the genetic network will be characterized by one of 
these steady states. When the signal generator starts to oscillate with an amplitude that doesn't 
drive the network far away from its steady state, the linear model is a good approximation. For 
large amplitudes, the nonlinear effects start to be important, and at some values of the generator's 
amplitude, the network will jump close to a different steady state. Such nonlinear behaviors can 
not be described by a linear model. Also, the parameters that describe the network are supposed to 
be constant in time. This approximation is valid if the changes in the network parameters are slow 
with respect to the changes produced by the oscillatory input signals. The input frequency should 
be chosen so that the system can be considered with constant coefficients for the elapsed time of 
measurements. Also, the period of oscillations must be less than the trend effects due to growth, 



13 



apoptosis etc. Beside biological effects that span large intervals of time, experimental artifacts, like 
medium evaporation, can superimpose a trend on the measured profile. The input period should 
be less than the time characteristics of these trends. This will impose a limit for the lower range of 
the input frequencies. The response to oscillations depends also on the time characteristics of the 
system under study. If the system has a high damping factor, the high frequencies will be strongly 
attenuated and the output signal is not measurable. With all these restrictions, the experimentalist 
still has the freedom to work in a frequency band, a freedom not present in the step stimulus. 

A different line of thought emerges when it comes to analyzing if the oscillatory method can be 
scaled to large networks. Experimentally, using high throughput measurements (microarray and 
proteomic tools) a large set of gene products can simultaneously be measured. The experimentalist 
is searching for a pathway that is controlled by a gene. Using oscillatory signals to stimulate the 
desired gene, the time variation of the downstream genes will contain in its spectrum the input 
frequency and so these genes will be detected. Moving the signal generator along the pathway, 
more and more local patches of the network will be uncovered. The global view on the network 
will consist of all these patches connected together. The theoretical framework for connecting a 
set of patches is unclear to us at present. Experimentally however, we verified that a source of 
oscillations propagates into a large genetic network, Specifically, a microarray experiment was 
conducted on mice entrained for two weeks on a 24 hours period of light-dark signals. The periodic 
input signal was not implemented at the level of gene promoter; it was an exterior periodic source 
of light that entrained the internal clock of the cell. After entrainment, and in complete darkness, 
the output signals (mRNA) were measured every 4 hours for 2 days using and Affymetrix platform. 
From about 6000 expressed genes in heart, about 500 showed a mRNA that oscillates with a 24 
hours period. Same results were reported in [23]. The next step is to implement the generator at 
the gene promoter level, and measure the spread of the input signal into the network. 

Given the advantages of a periodic stimulus presented above, we believe that the experimental 
implementation of a periodic generator at promoter level will prove fruitful in the study of genetic 
networks. 

Appendix 

The genetic network is described by a linear stochastic network pP, [Tl|, The network is 
driven using signal generators placed inside the promoters of a subset of genes that are part 
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of the network. For a gene we will denote by (D,r,p) the number of DNA, mRNA and pro- 
tein molecules respectively per cell. We consider r,p variables but D constant and we nor- 
malize it to D = 1. The state of a cell that contains n active genes is specified by: q = 
(Di, D2, D n , n, T2, r n ,pi,p2, ...,p n ). The genetic state is changing in time; for a short transi- 
tion time, dt, only one q~i changes its value and this new value can be either q~i + 1 or qi — 1. We 
consider in this paper a linear stochastic genetic network characterized by the following transition 
probabilities: T(q;q + lf,t) = YljLi AijQjdt, T(q;q — U;t) = ^2jLi^ijQjdt- Here q is the initial 
state and lj is a vector of length M with all elements except the one in the position i which is 
1. The time variation of the generators that drive the genes' expressions are encapsulated in the 
matrix A^ which governs the production of different molecules. The matrices A^ and Tij consist of 
four submatrices, corresponding to splitting the state q in two subgroups. One subgroup contains 
only the DNA states (D\, ■ ■ ■ , D n ) and the other subgroup contains the protein and mRNA states 
q = (ri,r 2: ...,r n ,p 1 ,p 2 ,...,Pn), 



A 





















, f = 






Gen 


A 







r 



(7) 



The generator submatrix Gen has a special form. It is a 2n x n matrix and locates the posi- 
tion of the generators in the genetic network: Genij = gi(t)5ij, i = 1 . . . 2n, j = 1 . . . n. Each 
gene promoter is driven by one generator gi(t),i = 1, ...,n, which will influence the mRNA 
production of gene i. The same mRNA production can be influenced by the protein concentra- 
tion, and this feedback effect is described by the elements of the 2n x 2n matrix A, (|7|). The 
structure of the matrix A is a consequence of the topology of the genetic network. The equa- 
tion for the probability P(q,t) of the network to be in the state q at time t is: dP (q,t)/dt = 
YhLi { E i~ ~ 1) Efeli A ik q k P (q, t) + YhU { E t ~ *) YlkLi fikQkP (q, t) , where the shift operators 
Ef are given by EfP(q,t) = P(qi, ± 1, ...Am)- 

We need the time evolution equations for mRNAs and proteins: fa =< qi > and =< qiqj > 
— < qi >< qj >, i, j = 1, . . . 2n. In matrix notation, for the column vector u and for the matrix X 
with elements given by Xij = — 5ijUi we obtain: 
d 



dt 
d_ 

~dt 



u = Hu + G , 

X = HX + XH T + Hdiag (u) + diag (u) H T + 2diag(T/j) . 

15 



(8) 
(9) 



Here H T is the transpose matrix of H = A — F and diag(n) has nonzero elements only on the 
principal diagonal: diag(fi)ij = Sijiii- Using the Laplace transform, the solution to (JHJ) is (1). 
The second equation, ©, is a matrix equation. To solve this equation we first transform it to 
an equation were the unknown is a column vector. The transformation needed is I h vec(X), 
where the column vector vec(X) contains the columns of the matrix X one on top of the next one, 
starting with the first column and ending with the last column. The vec mapping has the useful 
property that vec(HX) = (1 <g) H)vec(X), vec(XH) = (H T <8> l)vec(X), were 1 is the unit matrix 
and A <g) B is the tensor product of two matrices A and B. The column vector vec(diag (/j,)) can 
be expressed in terms of the column vector fx: vec(diag (//)) = Lfi, were L is a lift matrix from a 
space of dimension of \i to the square of this dimension: L = (Pi, . . . , P2n) 7 \ {Pk)ij = ^ik^jk- The 
solution to @ takes the form (2). 
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A Supporting Material 



A.l Derivation of the time evolution equations for the mean and fluctuation 
driven by signal generators 

The genetic state q = (D\, D2, D n , r±, T2, r n ,Pi,P2, ■■■■,Vn) is changing in time; let the state 
be qinitial at time t\ and q final a t a later time £2- The probability of transition from the initial 
to the final state is, in the most general case, a function of the initial state, the final state and 
the times of transitions: T(q~i n itiai',q~final',tiit2)- Following a common hypothesis, the transition 
probability is proportional with the transition time — 1\ if this is very short (£2 =t\+dt with dt 
an infinitesimal small quantity). The transition time being short only one g, changes its value and 
this new value can be either t/j + 1 or qi — 1. We consider in this paper a linear stochastic genetic 
network characterized by the following transition probabilities: T(q; q+li, ; t, t+dt) = Y^fL\ AijQjdt, 
T(q; q — If, t, t + dt) = X^=i ^ijQjdt. Here q is the initial state and lj is a vector of length M with 
all elements except the one in the position i which is 1. The equation for the probability of the 
network to be in the state q at time t, P(q,t), is then, |2j,: 

„ M M M M 

-P (q, t) = J2 (EC - 1) ~ A ^ P & *) + E i E t - 1) E ^QkP (q, t) , (10) 

i=l k=l i=l k=l 

where the shift operators are given by 

EfP(q,t)=P(qi,...,qi±l,...,q M ). (11) 

We need to obtain the time evolution equations for < q a > and < q a qp >, 

00 

< q a >=Y,q^ p (q^) > (12) 

00 

< q a q/3 >= E q^qpP (q, t) ■ (13) 

q=0 

The easiest way to follow the computations is to use the z-transform of a function, defined by: 

oo 

Z(P(q,t))= ...z M ^ M P(q,t) . (14) 

qi=0,...,q M =0 

The argument z of the z-transform will be displayed using the notation 
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F(z,t)=Z(P(q,t)). (15) 
The quantities of interest are related with the z-transform through 

F a = < q a > , (16) 
Fa/3 = < q a q/3 > -8 a {3 <q a > ■ 

where 5 a /3 is the Kronecker delta symbol which is if a ^ j3 and 1 if a = (3 and 

F a = ^ F (z,t) U= m =i...m, (17) 
Fa, = ^^(M)U=M=,..M • (18) 

The derivatives of the z-transform are not directly related to the covariance matrix: 

Va/3 =< qaq/3 > ~ < <7a >< Q{3 > ■ (19) 

However, the covariance matrix can be easily expressed in terms of the z-transform variables: 

z-a/3 = F a p - F a Fj3 + 5a/3F a . (20) 

The equation for F can be obtained by taking the z-transform of the master equation (jlUp using 
the following rules: 

Z(EfP(q,t)) = zr l Z(P(q,t))-zr l Z(P(q,t)\^ =0 ) , (21) 
Z(Erp(q,t)) = Zi Z(P(q,t)) , (22) 
Z(ftP(g,i)) = Zi d Zi Z(P(q,t)). (23) 

(24) 

If the degradation matrix T is diagonal pQ, then the probability P(q,t) \q i= o of the state with 
a missing molecular specie will not be part of the the equation for the z-transform. Indeed, the 
boundary term in the z-transform of EfTaqiP{q,t) will vanish for qi = 0. For a non-diagonal 
r matrix, we obtain the same equation if we work with natural boundary conditions, that is 
P(q,t) = if qi = for one i from the set 1...M. The majority of the genetic networks will not 
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obey the natural boundary conditions. However, the final results are the same for (i) a non-diagonal 
r matrix with natural boundary conditions and (ii) a diagonal T matrix with no restriction imposed 
on the boundary. For the sake of the symmetry of the computations, we will derive the results for 
a general T matrix and natural boundary conditions, and use a diagonal T matrix when we will 
study the behavior of a genetic network. 



The equation for the z-transform now reads: 

q M M q M M Q 

—F (z, t) = Y^ ( z i ~ 1) A *kZk-^F (z, t)+J2 - 1) H FikZk-Q-F (z, t) . (25) 

i=l k=l Zh t=l fc=l Zk 

Take the derivative of these equation with respect to z a 



i,k=l 

r ik l-Zi- 2 5 ia z k —F(z,t) + (zr 1 - 1) 5 ka —F(z,t) + (zr 1 - l) z k Q Qz F{z,t)\ . 

i,k=l 

Introducing Zi = 1, i = 1...M we obtain the equation for the time evolution of the mean values: 

d M 

-F a = Yi^ak ~ f ak )F k . (26) 
k=i 

For the second moments we continue to take derivatives of J2 



d 3 M 

dtdz dzp F ^' ^ = ^ A ik {5i a 8 kl3 d k F + 5 ia z k d kf3 F + 5 i(3 5 ka d k F + (z { - 1) 5 ka d kf3 F + 



i,fc=i 



$ipz k d ka F + (zj - 1) 5 k/3 d ka F + - 1) z k d kaf3 F) + 
^ f ik (2zi~ 3 5ii35 ia z k d k F - Zi~ 2 5 k f35 ia d k F - zC 2 b iOL z k d kf3 F - 

zr 2 5 il3 5 ka d k F + (zj -1 - 1) S ka d k/3 F - Zi~ 2 5 if3 z k d ka F + 
(zr 1 - l) 5 kf sd ka F + (zj -1 - l) z k d ka pF) , 



»,&=! 



d 3 



M M 



-F (z,t) \ Zi =l,i=l..M — A aP dpF + ^ A ak d k pF + ApadaF + ^ A pk d ka F + 
Q " fc=i fc=i 

M Af A/ 

^ 2 rp k d k F5 a p — F a pdpF - ^ F ak d k pF — rp a d a F — ^ rp k d ka F , 
fe=i fe=i fc=i 
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d M M 

~Tj ? ap = ^2 — r a kj F k p + ^2 (^A-pk ~ Fpkj F ka + AapFp + Ap a F a (27) 

k=l k=l 

M 



-TapFp — PpaFa + 2 5 a /3 ^ Ffi k F k . 



k=l 



This is the equation that we need. Later we will use it to reveal the action of the generators, 
that are hidden now in the coefficients Ai k . Before we deal with the generators, we will derive a 
general formula for the covariance matrix v a p to see how different it is from the one above. 

v a p =< q a qp > ~ < q a >< qp >= F a /3 ~ F a F(3 + S afS F a , (28) 

= - {jt F «) * - F 4t F ? + 6 °4 F " • (29) 

Now we insert the derivatives for F a and F a p 



d M M 

~^ Ua P = ^2 {^ak — r a kj F k p + ^2 {Apk — ^SfcJ F ka + AapFp + Ap a F a (30) 

k=l k=l 

M 



—FapFp — r/3 a F a + 2 5 a p Fp k F k + 

k=i 

M 

^2(-A ak F k Fp + f ak F k F p - Ap k F k F a + f pk F k F a ) + 

k=l 

M 

s a /3 ^2 yAak — r a k^j F k . 



k=l 

We want to get rid of the variables F a p and write everything in terms of v a p and < q a > . First 
we regroup the terms and then add and subtract the term 

M 

J2 (4** " F «k) hpFp (31) 
fc=l 
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to obtain 

d M M 

jfap = { A <*k ~ r ak ) (F kf3 - F k F p + 6 k(3 F p ) - (X* ~ F ak ) 8 kp Fp + (32) 

k=l k=l 

AI M 

+ ^ (As* - ffrk) ( F ka - F k F a + 5 ka F a ) - Y (Apk - f 0k ) 5 ka F a + 

k=l k=l 

M 

+A a pFp + Ap a F a - r a pFp — Fp a F a + 2 5 a p ^ r^ k F k + 

fc=i 

M 

+5aj3 ^2 (A~ak ~ F ak ) F k 



k=l 



d M M 

~^ y ^P = 5^ {j^ak — f ak j V k /3 + ^ \Aj3k — Ppkj V ka — (33) 
k=l k=l 

—AapFp + r a pFp — Ap a F a + r/3 a F a + 
+A a f3Fp + Ap a F a — r a pFj3 - rp a F a + 

AI M 

+<W ^ F ak F k + 5 a p ^ A ak F k 

k=l k=l 




A. 2 The Generators 

The generators constitute a submatrix of the matrix A : 

A=\ (35) 

Here Q/ g and a b are matrices with all elements zeros, where a, (3 = 1 . . . n, and a,b = n+1 . . . 3n. 
The matrix Q a p contains the generators and thus is a matrix with time dependant elements. The 
matrix A a b has constant elements which depend on the genetic network. From now on we make 
a distinction between Greek indices and Latin indices, so that we can rewrite the general time 
dependance equations in terms of generators. The Greek indices run along the DNA variables, 
whereas the Latin indices run through the mRNAs and proteins variables. 
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Lets specialize the equation (|27j) for Latin indices. We will split the summations using a generic 
Greek letter 7 and a generic Latin letter g. We consider the number of DNA be constant in time 
and normalized to the value 1. As a consequence 



Fjb =< Qjqb > -<V < % >=< qjQb > = < Mb >=< Qb >= Fb (36) 
In terms of Greek and Latin indices, the matrix T, looks like 



Q/ 3 a b 
Oa/3 r ab 



(37) 



so 



~faFah — 52 (A ay — ^07) F lb + 52 ( Aa 9 — F ag ) F gb + 52 (A bl — r &7 ) F ia + 52 (A bg — r bg ) F ga + 



+A ab F b + A ha F a - r ab F b - r ba F a + 25 ab (52 a 7 f 7 + 52 r bg F g j , 

V 7 9 I 



jFab = (%2 Ga ^ <Qb> + E Gb ^ <^>+52 ( Aa 9 - Fa 9) Fgb + 52 ( A bg ~ r bg) Fga + 

119 9 

+A ab F b + A ba F a — r ab F b — r ba F a + 2 5 ab 52 Fb g F g . 

9 

We have to eliminate the sum of the generators. We use for this the equation for the mean, 
taking care that for DNA variables, F 7 =< q 7 >= 1 



j t F a = (4ry " ^7) ^7 + 52 ( Aa 9 F 9 ~ r ag) F g (38) 

7 9 

= ^9^+52^9- Fag) F g . 
1 9 

We obtain then: 

jF ab = (j- f F a - 52 (Aag - r ag ) F g ^j F b + (j^Fj, - 52 (A bg - r bg ) F g ^j F a + 

+ 52 ( A ag - Fag) F gb + (A bg - r bg ) F ga + A ab F b + A ba F a ~ 
9 9 

— F a bFb — r ba F a + 2 5 ab 52 FbgFg , 

9 
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Jt Fab = (Jt Fa ) Fb + (Jt Fb ) Fa + ^ {Aa9 - Fa9) {Fab - FaFb) + s ^ - Fb ^ ^ - w + 

^ ' ^ ' 9 9 

+A ab F b + A ba F a - r ab F b - r ba F a + 2 6 ab ^2 r bgFg ■ 

9 

From the formula above, we see that a new variable appeared in a natural way: 

X ab = F ab - F a F b . (39) 
The time evolution of this new variable is given by the equation: 

j f X ab = ^2(H ag X gb + H bg X ga ) + A ab <q b > +A ba <q a > - (40) 

9 

-Tab <%> ~r ba < q a > +2 S ab ^2r bg <q g > , 

9 

with 

H ab = A ab - F ab , (41) 

or 



^-X ab = ^2(H ag X gb + H bg X ga ) + H ab <q b > +H ba < q a > +2 5 ab r bg < q g > . (42) 

9 9 



In what follows we will use a diagonal T matrix. For this case the equation simplifies to 



jX ab = ^2(H ag X gb + H bg X ga ) + A ab <q b > +A ba < q a > . (43) 

9 



The meaning of the matrix X can be found if we write it in terms of the covariance matrix u ab . 

Xab = F ab - F a F b = ~ $abF a = v ab - 5 a b < q a > • (44) 
Thus X measure the deviation of the stochastic process from a Poissonian process, 
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Vab = S ab < Qa> +X ab . (45) 

Now it is easy to write everything in terms of the reduced state q = (ri,r2, ...,r n ,pi,p2, ...,p n ). 
To do this we observe that q k = qk+m k = 1 . . . In. In other words, we subtract n from from each 
Latin index and keep the same notations for the variables. We use i = a — n, j = b — n, k = g — n. 

First, for the equation for the mean we simplify a relation deduced before, (|35|) 



-jj_Fi — ^2 (^«7 — A7) F-y + ^ {AikFk - r ik ) F k 

7 k 

= y~! + y~] (Ant - A*;) jpfc • 



(46) 



7 fe 

Note that in the sum Qi^ only one term is nonzero for i = 1 . . . n and all terms are zero for 
% = n = 1 . . . 2n. Indeed, each mRNA is controlled by only one generator: 



717 



$iygi(t) ■ 



(47) 



The above formulas tell also that only the mRNA is under the control of the generator, not the 
proteins neither the DNA. To simplify the notation we will write 



Gi{t) = gi(t), i = l...n, 
Gt(t) = 0, i = n + 1...2n 



(48) 
(49) 



The equation for the mean then simplifies to 



dt 



<Qi>=Y^ Hik <( lk> +Gi(t) . (50) 



A. 3 Solution to the Mean and Fluctuation equations 

The two equations from the previous section can now be written using a matrix notation: 



d_ 

d 



dt 



[i = H/i + G , 

X = HX + XH T + Hdiag (ji) + diag (p) H T + 2diag{Tu) , 



(51) 
(52) 
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where the column vector fi has the components fa =< q% >, and the matrix X is related with 
the covariance matrix v by: z/j,- = Sij < q\ > +X{j, i,j = l... 2n. Here H T is the transpose matrix 
of H = A — r and diag(fi) has nonzero elements only on the principal diagonal: diag(fi)ij = 5ijfa. 
We took care of the fact that X is a symmetric matrix, X T = X. 

The first equation in (|51|) has a column vector as an unknown, fa and is easy to solve it if we 
use the Laplace transform 

/•oo 

= / e - st fj,(t)dt . (53) 
J o 

The equation for the mean becomes 

s/i(s) - hq = Hfx(s) + G(s) , (54) 

with fiQ being the value of the mean number of molecules at time zero, when the generator was 
applied. Thus: 

fi(s) = (s-H)- 1 (G(s) + fi ) (55) 

The next goal is to solve for X. The second equation in (j51|) is a matrix equation. To find 
a solution for X we transform the matrix equation into a vector equation. The transformation 
needed is (jS] page 244): 



X i-> vec(X) , (56) 

where the column vector vec(X) contains the columns of the matrix X one on top of the next 
one, starting with the first column and ending with the last column. In index notations, the element 
Xij of the matrix X gets into the line i + m(j — 1) in vec(X) if X is an m x m matrix. 

The vec mapping has the useful property that 

vec(HX) = (1 ® H)vec(X) , (57) 
vec(XH) = (H T l)vec(X) , (58) 

were 1 is the unit matrix and A B is the tensor product of two matrices A and B. The matrix 
A (g> B is constructed by substituting each element a^- of the matrix A by the matrix aijB. 
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The matrix equation for X becomes 



d_ 

dt 



vec(X) = (H ® 1 + 1 <g> H)vec(X) + (# <g> 1 + 1 ® J?) vec(diag (pi)) + 2vec(diag(T u)) . (59) 



The column vector vec(diag (//)) can be expressed in terms of the column vector \x: 



vec(diag (fj,)) = Lu 



(60) 



were L is a lift matrix from a space of dimension of \x to the square of this dimension. The 
matrix L has the block structure 



L 



V P2n J 



(61) 



where 2n is the number of rows in fi ( n rows for mRNA and another n for proteins). The 
submatrices Pk, k = 1...2n are In x 2n square projection matrices, with all elements zero except 
one: 

(Pk)ab = dakhk ■ (62) 

As an example, for n = 1 we have 1 mRNA and 1 protein and the dimension of L is 4 x 2 



L 



With the same lift matrix L we can write 








(63) 



vec(diag(T fj,)) = LT^i 



Denote now the Laplace transform of vec(X) as V . We have, from 1591 



(64) 



sV(s) - V Q = (H®l + l® H)V(s) + ((H (g)l + l<g)iT)L + 2LT) n(s) , 



(65) 
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V(s) = (s- l®H-H®iy 1 ((l®H + H®l)L + 2Lr)(s-H)- 1 (G( y s)+no) 



(66) 



For a diagonal T 

V(s) = (s-l®H-H®l)- 1 (A®l + l®A)L(s-H)- 1 (G(s)+fi ) + 
+ («- 1 ®H-H® 1)~ X V , 

with no is the initial condition for the mean and Vo the initial condition for vec(X). 
From the above formula (|55|) we see that the mean values are expressed in terms of the generators 
through the mean transfer matrix: 



1 



(67) 
s — H 

The interesting form, (|66jh is the fluctuation transfer matrix, that passes the time variation of 
the input generators into the time variation of vec(X): 



s- 1 (8) H - H (8) 1 



[{1 ® H + H® l)L + 2Lr] 



s-H 



For a diagonal T matrix this simplifies to 



s - I® H - H ®1 
As an example, if H and V are 2 by 2 matrices 



[(1 ® A + A<& 1) L] 



s-H 



H = 


hu 




, r = 


5n 


512 






h,22 




_ 521 


522 



(69) 
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we get: 



s - 1 <g> if - if <g> 1 



-2 h n + s 
-h 2 i 
-h 2 i 




(1 ® H + H ® l)L + 2Lr 



-hi2 -hu 

-/in - h 2 2 + s -hu 

-/in - h-22 + s -h\2 

—h 2 i —h 2 i -2h,22 + s 

2/iii + 2t/n 2^i2 
/121 hi2 
h-21 hi t 2 

2 521 2/l 22 + 2 ^22 



(70) 



(71) 



((i ®H + H® l)L + 2Lr) 



s-H 



1 

A 



2/ins - 2/in/i 22 + 2 511s - 2511/122 + 2512/121 2/112/in +2/112511 + 2512s - 2512/in 

h-21 (s - h 2 2 + h 12 ) hi2 (/121 + s - /in) 

/l21 (s - /l22 + /ll2) /il2 (/l21 + S - /ill) 

2521S - 2521/122 + 2/121/122 + 2/121522 2521/112 + 2/l 22 s - 2/111/122 + 2522S - 2522/111 



A = s - h 2 2S - hns + hnh 2 2 - hi 2 h 



21 



(72) 



A. 4 An autoregulatory gene with a periodically driven cofactor. Response of 
the system to an arbitrary input 

One of the most fundamental regulatory motif in a genetic network is an autoregulatory gene 
through a negative feedback, [I] . We consider the case when the gene regulation is under the control 
of its own protein product and the protein activity is modulated by a cofactor. The equation for 
the mean is: 

"7r 



d 


' (r) ' 




lit 


. <P> . 







-h 

kp ~1p 





' (r) ' 




k + g(t) 






+ 






. (P) . 








(73) 
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where the state is q = (r,p). The cofactor, represented here by the term g(t), is driven by the light 
generator. The cofactor modulates the mRNA mean number (r) through an additive coupling. For 
this case we use suggestive notations Xn = X rr , X\i = X rp etc.. The Laplace transform of the 
quantities of interest, ((r(i)), (p(t)) , X rr (t), X rp (x), X pp (t), g(t)) are denoted by the same letter but 
the argument being the complex frequency s instead of the time t, like in 

poo 

(r)(s) = / (r(t))e- st dt . (74) 
J o 

The values of the mean number of molecules and their fluctuation, will depend on the inter- 
nal parameters j r ,^ p ,h,k p ,kQ as well as on the external parameters of the generator g(t). Two 
important natural parameters of the system play a significant role: 

= hk p + 7 r 7 p , (75) 
= 7r + 7 P • (76) 

The mean number of molecules are connected to the generator through: 



" (r) (s) ' 


1 


S + lp 


-h 




' 9(s)~ 


_ (p) w . 


ago 










k p 


S+Jr 








(77) 



with 

A (s) = s 2 + sujx + ujq 2 . (78) 

The deviation from a Poisson process measured by the variable X is under the generator influ- 
ence also: 



X rr (s) 




-2h(s + 2j p ) k p (s + 7 P - h) 


2h 2 (s + 2 7p ) (s + kp + jr) 


X rp (s) 


1 


k p (s + 2j p ) (s + 2 7 r ) (s + 7 P - h) 


-h (s + 27 P ) (s + 2 7 r ) (s + k p + 7 r ) 


X pr (s) 


A f (s) 


kp (s + 2 7 p ) (s + 2 j r ) (s + 7 P - h) 


-/i (s + 27 p ) (s + 2 7 r ) (s + k p + j r ) 


X pp (s) 




2k p 2 (s + 2 7r )(s + 7 P - h) 


—2h(s + 2j r )kp(s + kp + j r ) 



(79) 

with 

A f (s) = (s + wi) (s 2 + swi + uj q 2 ) (s 2 + 2 swi + 4w 2 ) (80) 
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A. 5 The step and the periodic stimuli 



There are two cases of interest to us, a step stimulus and a periodic one. 
For a step stimulus: 

ow = £. 



(81) 



We consider that the system is in a steady state before we apply the step stimulus. The steady 
state is governed by the translation rate ko. For a stable system ( Re(Ai.2) > 0), the mean number 
of molecules decay exponentially to zero. 



A1A2 



A1A2 Ai (A2 — Ai) 



(Ai-7p) Ge -Xit + (A2~7p) Cc -\it 
A2 (Ai — A2) 



ki 



G + 



kr; 



Ge~ Xlt -f 



km 



A1A2 u A2A1 X\ (Ai — A2) A2 (A2 — Ai) 

where Ai and A2 are the eigenvalues of H: A (s) = (s + Ai) (s + A2) , 



Ge 



(82) 
(83) 



Ai = 1/2 wi - 1/2 v^i 2 -4w 2 
A 2 = 1/2 wi + 1/2 ^wi 2 -4w 2 



(84) 
(85) 



For fluctuations we get also exponentially decaying responses to the step stimulus: 





(/-) 


= X rr fi 




-u>lt 


+ X rr) \ i e 


X\t 1 v — 
-t-s\ rr \ 2 z. 


Xot 1 \r „ — 2Ait 1 v — 2 A 2 i 
+ A rr, 2Ai e + A rr,2A 2 e j 




(t) 


= X rp fi 


-\- X r p W1 e 


-W\t 


+X rP) \ 1 e~ 


~\~X r p\ 2 & 


" A2i +^r P ,2Aie- 2Ali +X rPi2 A 2 e" 2A2i 




(/-) 


= Xppfl 


~^~Xpp t ( J j 1 e 


-U\t 


+X pp ,\ 1 e 


-Ait , y 


— X^t 1 v" Ait , \- „—2\'zt 
+ App,2Ai6 +A pPi 2A 2 e 
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The coefficients from the above formulas are collected in the following matrices 



X rr fi 




X r r,u>\ 




Xrr,Xi 




X rr ,x 2 




X rr ,2X! 




X rr ,2 A 2 





fefcpC7(fe-7p)7 p fefcpfco(fe-7p)7 P 



A 2 2 Ai 2 u;i 



+ 



A 2 2 Ai 2 o;i 



fcfcpG(-2 7p+ui)(-7 P +fo+ui) 

(-2Ai+o;i)(-A 2 +wi)(a;i-2A 2 )(-Ai+LJi)a;i 

_ 2 hk p G(- Ai +2 7p)(Ai+fe-7p) 
(Ai-2A 2 )(-A 2 +Ai)Ai 2 (-Ai+^i) 

_fj fefcpG(2 7p-A 2 )(A 2 +fe^7p) 
A 2 2 (-A 2 +Ai)(-A 2 +cji)(2Ai-A 2 ) 

fefcpG(-Ai+7p)(2Ai+fe-7 P ) 
Ai 2 (-A 2 +Ai)(2Ai-A 2 )(-2Ai+^i) 

fefcpG(7p— A 2 )(2 A 2 +fe~7p) 
A 2 2 (-A 2 +Ai)(Ai-2A 2 )(c;i-2A 2 ) 
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X 



rp,0 



X 



rp, uii 



-"^Vp,A 2 
X rp ,2 Ai 
X r p,2 A 2 



fcp(fco+G)7 r 7p(fe-7 P ) 



A 2 2 Ai 2 o;i 



fcpG(^i-2 7 r )(-2 7p+^i)(-7 P +fe+^i) 
wi(^i-2 A 2 )(-A 2 +o; 1 )(-Ai+wi)(-2 XT+^l) 

fcpG(-Ai+2 7r )(-Ai+2 7p)(Ai+/i-7p) 
Ai a (-A 2 +A 1 )(-Ai+c; 1 )(Ai-2A2) 

fc P G(-A2+2 7 r )(2 7p-A 2 )(A2+/i-7p) 
(-A 2 +o; 1 )(2Ai-A 2 )(-A 2 +Ai)A 2 2 

fcpG(-Ai+7 r )(-Ai+7p)(2Ai+/t-7 P ) 
A 1 2 (-A 2 +A 1 )(2A 1 -A 2 )(-2A 1 +^ 1 ) 

fcpG(-A 2 +7r)(7p-^2)(2 A 2 +fe-7p) 
(-A 2 +A 1 )(A 1 -2A 2 )(w 1 -2A 2 )A 2 2 







Xpp^u) x 




X ppM 








Xpp,2 Ai 




_ -Xpp,2A 2 
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fcp 2 (fcQ+G)7r(fe-7p) 

A 2 2 Ai 2 wi 

fcp 2 G(-7p+fe+(^i)(a;i-2 7 r ) 

wi(wi-2Aa)(-Aa+a;i)(-Ai+wi)(-2Ai-|-a;i) 

2 fc P 2 G(Ai+fo-7p)(-Ai+2 7 r ) 
Ai 2 (-A 2 +Ai)(Ai-2A 2 )(-Ai+lji) 

2 fcp 2 G(A 2 +fe-7p)(-A 2 +2 7 r ) 
(-A 2 +Ai)(2A 1 -A 2 )(-A 2 +o;i)A 2 2 

fcp 2 G(2 Ai+fe-7p)(-Ai+7 r ) 
A 1 2 (2A 1 -A 2 )(-A 2 +A 1 )(-2Ai+oji) 

fcp 2 G(2 A 2 +/i-7p)(-A 2 + 7r ) 
A 2 2 (Ai-2A 2 )(-A 2 +Ai)(oji-2A 2 ) 

For the periodic case with an input frequency to and amplitude a, g(t) = ko + acos(ut), ( fco is 
a baseline not controlled by the exterior light input) 

k 



9{s) 



+ 



s 2 + w 2 



(90) 



We keep only the stationary solutions in the response ( in practice we wait for the transients to 
become small enough) 



<r(t)> 

<P(*)) 
X rr (t) 

X rp (t) 

Xpr (t) 



Ro + Rxe^ 1 + R\e- iu}t , 
P + Pie iujt + P^e~ iujt , 
X r fi + X r%i \e lult + X*ie lw , 
X rp fi + X rp% i e* w ' + X* pl e luJ t 
Xp : o + Xp t ie lult + X*ie * w< . 



(91) 
(92) 
(93) 
(94) 
(95) 



The star * means complex conjugation. In terms of the parameters that constitutes the autoregu- 
latory system we have: 
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Ro = lp \ 7 , (96) 



Irlp + hk 
kpko 



p 



7 r 7 p + /ifcp 



a (7 P + iu) 
k„ a 



Ri = 1/2 » ( 98 ) 

wo — w z + zw w_; 



p i = 1/2 2 , . , (99) 

— W z + Wo + 100 LOl 

X r0 = Mft-7p)Myfe (100) 
wo Wl 
k Qlrlp k p ( 7p - h) 

wiw 4 

X Pi0 = ^(7p-^)7r (1Q2) 

W 4 W1 

x i = -za (-ij p + u + ih) (-w + 2 i-y p ) hk p 

(— w 2 + ujq 2 + &w wi) (— w 2 + 2 zw wi + 4 wo 2 ) (—a; + z'wi) ' 
_ 1 l0 ak p (w - i7 p + jfe) (w - 2 i 7r ) (w - 2 j 7p ) 

(W z — Wo — VjJUJ\) (OO z — 4 Wo - izWWi) (W - IU)\ ) 

x i _ ia (-Hp + uj + ih)(uj-2 ijr) k p 2 

p ' (w 2 — wo 2 — iw wi) (w 2 — 2 iw wi — 4 wo 2 ) (w — iuj\) 

The impact of the natural (internal) frequencies wo and wi on the protein and mRN levels and 
fluctuation can be read out from the absolute values of the denominators of the mean and X: 

A = \det{iuj - H)\ 2 = w 2 ((w 2 -w 2 ) 2 + w 2 W! 2 ) , (106) 

2.2/2, ,2\2 / / . 2 A . 2\2 . . 9 . 2 



A/ = |(iet(«w-l®i/- J H"®l)| 2 = 4wiWo 2 (wi +w 2 ) (^(w 2 -4w 2 ) +4w 2 wi J. (107) 

We observe that wo is a resonance for the mean and X), whereas 2wo is only for X. 

Beside the ratios expressed by formulas (5) and (6) from the main paper, we can form different 
combinations between the periodic response variables that become useful for estimating the order 
of magnitude of the coefficients k, /i,7 r ,7 P (we consider the case when no experimental noise is 
present): 



(108) 
(109) 
(110) 



l^ll 2 


1 2 
k*" + 


7p 2 


l^ll 2 
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h. p 


X r ,o 


hip 




Xpfi 


kplr 




Rl,w=Q = 


1 lp 

^tt a . 

2 w 2 
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Prom the first relation we can estimate k p and 7 P . From the second one we estimate the ratio 
The third equation gives the variation of mRNA amplitude with the input amplitude a for small 
u . From this relation we estimate ojq . The mRNA degradation coefficient j r can now be obtained 
from 

7r = ^o/(7p + — k p ) . (Ill) 

7r 

Now we have h from h/j r . The last parameter ko comes from 

Ro = ^- (112) 
There are other interesting ratios worth to be written down: 

1 k-n 



Pi,u=o = ^^«, (H3) 



|A rjl | 2 V^ 2 + 4 7r 2 ) 

P q k p ' 



(114) 
(115) 



These relations can be used to further verify the validity of the model, once we estimated the 
parameters. 

A. 6 Fluctuation resonance 

We want to find a driving frequency for which the fluctuations dominates over the mean values. 
For such a frequency the system will be in a pure fluctuation resonance. In such a situation the 
molecular noise can drive the cell out of its equilibrium state, which can have dramatic consequence 
on the cell fate. At the fluctuation resonance frequency, the deviation from a Poissonian process, 
measured by the quantity X, should be very high. To measure this deviation we consider the ratio 
of the fluctuation amplitude | X p i \ over the mean amplitude | P\ \ (an analog of the Fano factor 
in frequency domain): 

^J W , (^->>>W> \"\ (116) 



\ ((w 2 -4u;o 2 ) 2 + 4u;V) {uj 2 + uj^) 
For systems for which ujq 3> u>\ we can se a resonance for fluctuations but not for the mean 
values at the input frequency to = 2w - A plot of this ratio is presented in Fig.3. We notice that the 
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width and the hight of the resonance are inverse proportional. The parameters for which we see the 
resonance in Fig. 3 doesn't belong to case we studied for step stimulus (Ai^ are real numbers there 
and complex here). The response to the step stimulus for systems that can enter into fluctuation 
resonance is a superposition of damped oscillations. Even in this situation the transients are gone 
after few periods. 

A. 7 The Genetic Network Spectral Function 

The time response (mean and fluctuation) of the autoregulatory system to a step stimulus can be 
expressed in general as a sum of 6 terms 

fexpit) — S eX pfi ~\~ S eX p^\B + S e %p t 2e ^ 2 + S exp ^e ^ 3 ~\~ S exp ^e ^ 4 + S exPi §e ^° • (117) 

Only three of these terms are present in the mean. For the purpose of the following analysis, we will 
consider only the case when all rj's are positive, which is equivalent with uj\ > 2u>o. The asymptotic 
response, of the same autoregulatory system, to a periodic stimulus has the form 

fper (t) = S per fi + Sper,l e + ^per,l e ; (H8) 

for both the mean and the fluctuation. The parameters of the system k p ,h,~f r ,~f p are hidden in 
the coefficients S exp ^ or S per j, i = 0, . . . , 5, j = 0, 1. For more complex genetic network, the time 
evolution of the measured quantity f(t) can be expressed as 

f(t)= / S{x)K{xt)dx . (119) 

Here S(x) is the spectral function that contains the information about the genetic network and 
K(xt) is the kernel that depends only on the type of the stimulus ( i.e. on the experimental design). 
Indeed, for an autoregulatory network, using the Dirac's (^-function, we have 

fex P {t) = / S exp (x)e- Xt dx , (120) 

J a 

6 

Sexp (x) = S expji 5 (X - T]i) , (121) 

i=l 

K exp {xt) = e~ xt . (122) 
The values a and b are chosen such that the spectrum S exp (x) is zero outside the interval [a, b]. 
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For the periodic stimulus, we have a similar representation for the spectral function S(x) but 
the kernel is different 



fper (i) — / Sp er (x) e dx , 



(123) 



J-u 

Sper (x) = Sp er fiS (x) + S per ^5 (x — Lj) + S* per ^5 (x + u) 



(124) 



per \X 



(125) 



with Q > uj. 



The topology of the genetic network is reflected in the spectral function S(x). Given a set of 
measured data, first we have to recover the spectral function of the network and then from it the 
parameters of the network. If we lack a good model for the topology of the genetic network we 
cannot find the parameters of the network, but we can recover the spectral function S(x) from 
the data (the kernel K(xt) does not depend on the network). Thus different genetic networks 
can be compared using their spectral functions. However, the spectral function depends on the 
experimental design. We proved for the autoregulatory system that the spectral function S per is 
much simpler than S exp . We want to show that there is even a deeper difference between these two 
experimental designs. Namely, in the presence of experimental noise, it is much easier to recover 
S per from the experimental data than S exp . This phenomena appeared in other branches of science 
and in many different forms. To adapt it to biology, we noticed that a legitimate question from a 
molecular biologist is: instead of creating new assays to measure S per why is not enough to increase 
the number of replicates to obtain an accurate S exp ? We will prove that the number of replicates 
for S exp growth exponentially with the accuracy. In what follows we collect and use for our specific 
problem, results form 

In laboratory measurements, we don't have f(t) for all values of t. Rather, we have samples of 
it at discrete time points. For the periodic stimulus, we measure /(rn), where n = 0, 1, . . . , N. As a 
working example, consider the samples of the mean of the mRNA, r(n) = {r(rn)),n = ... N — 1, . 
The unknown spectrum S per {u)) and r(n) are related through the equation: 



There are three parameters in the problem: r,£l,N. The sampling parameter r must be such 




(126) 



3G 



that the input frequency uii n can be detected in the output data, that is r < n/ujin. The frequency 
£1 should be greater than the input frequency Ui n . There is no condition on the number of points 
N. Because we have a finite number N of measured data points, the spectrum S per (uj) can only be 
approximated as a weighted sum of N functions 3>fc(u;) ( see Appendix 1 and [H]) 

N-l 

S per {co) = s k ^ k (u). (127) 

fc=0 

The functions 3>fc(a;), k = • • • N — 1 come from a eigenvalue problem for an N x N matrix (see 
Appendix 1 at the end of this Supporting Material). Now, the experimental noise will alter the 
coefficients so the recovered spectrum will be: 



N-l , s 

S per (co) = + , (128) 

fc=o V PkJ 

where the e k are the noise coefficients. The numbers k = 1 • • • N — 1, come from the same 



eigenvalue problem as before and they depend only on the parameters r, N and not on the noise 
coefficients e^. Due to noise, we cannot use all N terms in (|128|) . but only the first J p , for which 

^<-, k = l,---,J p . (129) 

The right hand side of ()129|) is the Signal to Noise Ration (SNR) and for simplicity we will 
consider that is independent of the index k. The numbers decrease as k increase and so the 
condition for the cutoff J p is simple 

J_ < SNR < —L- . (130) 

PJ P P(J P +i) 

The exponential case can be developed parallel to the periodic case, 5 . The problem now reads 

like 

r (n) = / e~ p - x S exp (A) dX . (131) 



a 



Unlike for the periodic case, here a geometric sampling is optimum [S] 

Pn = -A n , n= 1...N . (132) 
a 

The limits a and b are chosen so that the spectrum is nonzero only inside [a, b]. For the periodic 
case we know the input frequency so we don't have to guess an interval [a, b] as we have to do 
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for the exponential case. Only the ratio 7 = b/a is important, as we see by changing the variable 
A = ax 

r (n) = J e~ apnX S exp {ax) adx . (133) 

Similar to the periodic case, solving an eigenvalue problem we can find an N-dimensional approxi- 
mation to the spectrum. Because of the experimental noise we can use only J e degrees of freedom, 
not N : 

S exp (A) = J2(s k + * fc (A) . (134) 

Here the terms e& are due to random experimental errors. Again, the functions \I/fc(A) and the 
numbers a k come from an eigenvalue problem ( different from the periodic one) and they don't 
depend on the noise but only on the parameters a, q, A, 7, N (actually, the numbers a& do not 
depend on the parameter a, only ^(A) does.) The cutoff J e is noise dependent and is given by 

— < SNR < . (135) 

a J e a (J e +i) 

The cutoffs J p and J e are of prime importance because they measure the number of degrees 
of freedom in the recovered spectrum. Desirable is that both cutoffs be as close as possible to 
the number of measurements, N, which is the case when the Signal to Noise Ratio (SNR) is 
high. Although the equations ()13Uj) and ()135|) look formally similar, they give completely different 
solutions to the cutoffs. This is a consequence of the different rate at which the numbers a k and 
(3k decrease to zero which we will study in the next section. 

A. 8 The number of replicates 

The SNR dictates how many spectral components are reliable and can be use to recover the spec- 
trum. We can imagine that by using replicates we can improve the SNR and so the two cases will 
come close to each other. This is not true; actually we need an unrealistic number of replicates to 
keep even few components for the exponential case. Indeed, with the help of r replicates, the SNR 
increase to 

SNRy/r , (136) 
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and the equations for the number of components J e , J p to enter into the recovered spectrum are 

— < SNR^ < — , (137) 

«J e «Je + l 

< SNRy/? < — L_ . (138) 
PJ P PJp+i 

The plots , Fig 7, of the number of replicates r as a function the number of spectral components J e 
or J p reveal that using a periodic stimulus we can use many more spectral components to recover 
the spectrum. The number of replicate growth very fast in the exponential case (for SNR = 10 
we need 269 replicates for 4 spectral components), whereas in the periodic case, the number of 
replicates stays low for many spectral components ( only for the 17th component it raises to 14, 
with SNR = 10). 

The source of such a discrepancy is that the eigenvalues a k tend fast to zero as 

« fc 2 = (139) 

COSh (7T&J 

where tends to infinity like a polynomial of degree at least one in k (there is no analytical formula 
for £ k ). For the plotted example, 7 = 5,q = 1/20, A = 60 1 / 20 
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10- 3 , 


(142) 


03 


= 1.74 • 


10- 5 , 


(143) 


ai9 


= 9.94 ■ 


IO" 29 . 


(144) 



For the periodic case the situation is much better. Here the numbers fit depends only on the 
product tVL and so is customary to introduce the parameter w through 2ttw = tQ. Then, for 
w = 1/3 for example, we get 



Ih 

02 

Ih 



0.99 , 
0.99 , 
0.99 , 
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(145) 
(146) 
(147) 



04 = 0.99 , (148) 
(3 20 = 0.000084 , (149) 

There is no an general analytical formula for /3k but it was proven that the first 2Nw beta numbers 
are close to 1 with the rest of them decreasing fast to zero. The fact that the majority of the 
eigenvalues for the periodic case are 1 whereas the eigenvalues for the exponential case decrease 
fast to zero is the source of the difference between the two cases. 




Signal to Noise Ratio 



Figure 5: The Threshold function of SNR 

Another interesting question is related to the resolution of the different exponentially decaying 
signals present in the output signal. For the periodic case we do not address this question, because 
the output signal has the same frequency as the input periodic signal (after the transients are gone) . 
However, for the response to a step stimulus, the transients contain the information. To obtain 
this information we have to resolve the transient components. The resolution power depends on the 
Signal to noise Ratio (SNR). For example, to be able to resolve the decay rates Ai and A2 when 
they are real positive numbers, we need to have 

— > Threshold(SNR) . (150) 

UJQ 

The Threshold as a function of SNR is plotted in the figure. Notice that we work with real Ai 2 
so uj\ > 2too for all SNR. 
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Appendix 1. The eigenvalue problem for the periodic case 

Recall that the measured quantities r(n) for n = . . . N — 1, can be expressed as 



r(n)= / e inTUJ S per {uj)cLu . (151) 
J-n 

From the N data points we can find a N-dimensional approximation to the spectrum S per (oj) solving 
the following singular value problem: find Vk(u)),Vk(n) and that satisfy 

CV k = \ k vk, (152) 
C*v k = X k V k , (153) 

where the operator C and its conjugate C* are 

(Cf)(n) = [ e in ™f(cj)ck;, (154) 
N-l 

{C*g){u) = Y, e ~ inTUJ 9(n) . (155) 

n=0 

The set Vk(u) form an orthonormal basis in L 2 (— 0,0) and Vf.(n) an orthonormal basis in the 
euclidian space E N . In the basis, the N-dimensional approximation to the spectrum reads like 

N-l 

sSLr M = E y v *M ' (156) 

k=Q k 

where the coefficients are obtained from the decomposition of the measured data r(n) 

N-l 



(n) = E r ^ v k (n) • (157) 



k=0 



The solution to the singular problem (|152j) can be reduced to the eigenvalue problem for the operator 
££*: find the eigenfunctions Vk{n) and the eigenvectors \\ from 



v— I sin (r (m — n)) , . ~ , N , 
E ? r^fcM = n , 158 



7r (m — n) 

m=0 ' 



where 



Pk = yi Xk - (159) 
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In this way, the solution to our problem is reduced to the diagonalization of an N x N matrix. 
This is the famous problem, 0. We have two independent parameters r and f2. The eigenvalues 
of the problem (|158|) depends on w, defined as 2irw = rf2. The first 2Nw eigenvalues are close 
to 1 with the rest of them close to zero. As a consequence, from H156|) we see that we can keep 
only the first 2Nw terms, because the rest of them are highly amplified by the small values of the 
eigenvalues which is dramatic when the values r k are corrupted by noise. We want than 2Nw to be 
close to N which case w = 1/2 and = tt/t. This situation corresponds to a sampling parameter 
r tuned for recovering the spectrum up to the frequency In general case when f2 > tt/t. The 
recovered spectrum, when noise is present will be than 

JV-l 



k=0 



^>)=£^^ fc M. (160) 



To connect with the notations from Section 10, denote &k( u ) = \l T/^V k (uj) and s k = r k /(3 k . 
Appendix 2. The eigenvalue problem for the step stimulus 

The problem for the exponential decay responses was solved in [5]. The unknown spectrum 
S exp and the measured data r(n) are connected through the equation 

r (n) = J\- apnX S exp (ax) adx , (161) 

with 7 = b/a and p n = (q/a)A n ,n = 1, . . . , N. Like for the periodic case, an N-dimensional 
approximation to the spectrum can be found from the solutions of two coupled equations: 

K,U k = a k u k , (162) 

K*u k = a k U k , (163) 

where 

(£/)(«) = J\- ap " x f(x)dx, (164) 
N 

{K*g){x) = Y J W n g(n)e- a ^ x , (165) 

n=l 

with the weights given by w n = p n ln(A), see [3]. The unknowns are the functions U k {x) that form 
an orthonormal basis in L 2 (l,j) and the functions u k (n) that form a basis in the euclidian space 
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M. N endowed with the scalar product 

N 

(g,h) = ^2w n g(n)h(n) . (166) 

n=l 

The N-dimensional approximation to the unknown spectrum can now be written as a decomposition 
in Uk basis as 

1 N 

S^ p {\) = -Y, r ^U k {\/a) , (167) 
fc=l 

with the components obtained from decomposing the measured data r n in the basis Uk 

N 

rk = ^2 w n r (n) u k (n) . (168) 

71=1 

Similar to the periodic case, the eigenvalue problem to be solved now is 



J- . /in in — p -a{Pn+Pm) _ f>-b(p n +Pm) 

> , 1 Uk{m) = a k u k (n) (169) 

with Uk(n) = ^Jw n Uk(n). The matrix that is diagonalized in the problem (|169|) is a symmetrized 
version of OC* and so there is a scaling difference between Uk and Ufc. The eigenvalues otk tend 
fast to zero as 

afc2 = J . n . ( 17 °) 

cosh(vr^) 

where tends to infinity like a polynomial of degree at least one. The recovered spectrum is 

1 N 4- 

where the terms are due to random experimental errors. 

In Section 10 we write the spectrum in terms of ^fe(A) = (l/a)Uk(X/a) and = rk/otk- 
Appendix 3. The eigenvalue problem for continuous measurements 
We discussed the spectrum recovery from a finite number of data, which is the case of laboratory 
measurements. However, it is instructive to inspect the case when we know f(t) from (|119|) for all t 
and in the limit for which a = 0, b = oo and Q, = oo. This problem was studied in jjj. As a bonus, 
we get an expression for the resolution of the exponential spectrum and a direct understanding of 
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the difference in the two eigenvalue problems presented in Appendix 3 and 4. The solution for the 
exponential case is in terms of the eigenvalues and eigenfunctions of the kernel K{^t) 

poo 

/ K(r ] t)<S> n (r ] )di 1 = Z n <£> n (t) . (172) 
J o 

The eigenfunctions form an orthogonal basis and both the measured data f exp (t) and the unknown 
S exp (ii) can be decomposed as: 

oo 

fex P (t) = J2fe*P*®k(t) , (173) 
k=l 

oo 

Sexp{ri) = ^2S exPtk ^ k {r]) . (174) 

k=l 

Now we have the relation between the spectrum and the measured data 

oo „ 

S(V) = ' ( 175 ) 

k=l ^ k 

with Hfc arranged in decreasing order Hi > H2 > .... We see from this expression that if 
decrease to zero and the components of the measured data f k are corrupted by noise, than the 
components with large k cannot be used to recover S exp (r]). The function thus recovered S exp (rj) 
has information just from the first components f eX p,k- Only if the eigenvalues don't decrease to 
zero we can use all the terms in the decomposition. 

For the exponential decay problem (|12U|) the eigenvalues form a continuous spectrum ( k is a 
positive real number) 

^ k ^ cosh (-7T k) ^ 
For the periodic solution ()123|) with Q = 00 (Fourier transform) the spectrum is discrete ( k = 
0, 1, ... ,00) 

E k = -i k V2^ (177) 

It is obvious the difference between the exponential decay situation ( step stimulus) and the 
periodic response. In the former case the eigenvalues tend fast to zero whereas in the later case 
they never approach zero ( they have a constant modulus one.) 

We aim now to find the resolution limit for resolving the exponential decay problem 0. For 
a given signal to noise ration SNR we want to find the minimum ratio of the exponential decay 
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rates rji/rji+i that can be resolved. Going a little deeper into the solution of the exponential decay 
response, j7] we find that the decomposition of the spectrum g(fi) is 

Sexp (r?) = / a h + E k + $> k + (77) dk + / a k ~ E k ~ (r?) dk (178) 
io Jo 

where the eigenfunctions are 

<f> k +(ri) = cos ( Hn (t?) - ^ ) (179) 

yfc vr \ ^ / 

$ fe -(r?) = --LsinjHnO?)-^ (180) 
VK vr \ ^ / 

with the angle expressed in terms of the Gamma function 

9 = angle (T (1/2 + ik)) . (181) 

Due to noise, we can recover the components up to a maximum fco, so we have all the components 
with k < ko. For this reason, we only can resolve points on the axis r\ that are separated at a 
distance larger than the distance between two zeros of $fc . Due to the presence of ln{r\) in the 
argumet of the trigonometric function, the zeros are 

1/2 e + m-K 

Mm = e fc o (182) 
To conclude, two decay rates rj a and can be recovered from the measured data if 

!h > jh^ = e (^) (183) 

Vb Vm+l 

The value ko that is the index for the maximum eigenvalue recoverable from noise is given by 
comparing the signal to noise ratio with the eigenvalue 

cosh (71" /cq) 



7T 



SNR 2 (184) 



Applying p83|) to the example we work with 1)84(1 we obtain the condition 



— > 2 cosh (185) 
loq \2k 
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